DEanalysis_kfish_osmotic
Lisa Johnson
2019-05-09
1 Introduction
A mixed model with random effects was chosen for this multifactor experiment, and analyzed using the limma package in R. This package implements a linear modeling approach and empiracal Bayes statistics. The limma package with the voom method estimates the mean-variance relationship of the log-counts, generates a precision weight for each observation and enters these into the limma empirical Bayes analysis pipeline.
In this model, clade (levels = Clade1, Clade2, Clade3), physiology (levels = Marine, Freshwater), and experimental condition (levels = 15ppt, 0.2ppt) are fixed effects while species (levels = 14) is considered a random effect.
The raw counts file, generated with NumReads from the salmon (version 0.12.0) quantification tool and summarized with the tximport Bioconductor package (version 1.10.1) in R, can be downloaded from an osf repository with this link, then imported into the R framework.
Samples from species with low numbers of replicates were dropped from the raw counts table (F. zebrinus, F. nottii, F. sciadicus). The raw counts table has the following dimensions (genes x samples).
#dim(counts_design)
#[1] 30471 130
# -----------------------
# Format design and counts matrix
# Drop columns with no data
# -----------------------
design <- counts_design[counts_design$Ensembl == 'Empty',]
#design$type <- c("species","native_salinity","clade","group","condition")
drops <- c("X","Ensembl",
"F_zebrinus_BW_1.quant","F_zebrinus_BW_2.quant",
"F_zebrinus_FW_1.quant","F_zebrinus_FW_2.quant",
"F_notti_FW_1.quant","F_notti_FW_2.quant",
"F_sciadicus_BW_1.quant","F_sciadicus_FW_1.quant","F_sciadicus_FW_2.quant")
transfer_drops <- c("F_sciadicus_transfer_1.quant","F_rathbuni_transfer_1.quant","F_rathbuni_transfer_2.quant","F_rathbuni_transfer_3.quant",
"F_grandis_transfer_1.quant","F_grandis_transfer_2.quant","F_grandis_transfer_3.quant",
"F_notatus_transfer_1.quant","F_notatus_transfer_2.quant","F_notatus_transfer_3.quant",
"F_parvapinis_transfer_1.quant","F_parvapinis_transfer_2.quant",
"L_goodei_transfer_1.quant","L_goodei_transfer_2.quant","L_goodei_transfer_3.quant",
"F_olivaceous_transfer_1.quant","F_olivaceous_transfer_2.quant",
"L_parva_transfer_1.quant","L_parva_transfer_2.quant","L_parva_transfer_3.quant",
"F_heteroclitusMDPP_transfer_1.quant","F_heteroclitusMDPP_transfer_2.quant","F_heteroclitusMDPP_transfer_3.quant",
"F_similis_transfer_1.quant","F_similis_transfer_2.quant","F_similis_transfer_3.quant",
"F_diaphanus_transfer_1.quant","F_diaphanus_transfer_2.quant",
"F_chrysotus_transfer_1.quant","F_chrysotus_transfer_2.quant",
"A_xenica_transfer_1.quant","A_xenica_transfer_2.quant","A_xenica_transfer_3.quant" ,
"F_catanatus_transfer_1.quant","F_catanatus_transfer_2.quant",
"F_heteroclitusMDPL_transfer_1.quant","F_heteroclitusMDPL_transfer_2.quant","F_heteroclitusMDPL_transfer_3.quant")
counts<-counts_design[!counts_design$Ensembl == 'Empty',]
rownames(counts)<-counts$Ensembl
design <- design[ , !(names(design) %in% drops)]
counts <- counts[ , !(names(counts) %in% drops)]
design <- design[ , !(names(design) %in% transfer_drops)]
counts <- counts[ , !(names(counts) %in% transfer_drops)]
#dim(design)
#[1] 5 81
dim(counts)## [1] 30466 81gene.names<-rownames(counts)
design[] <- lapply( design, factor)2 Sample Design Matrix
A matrix was generated using the following model:
~0 + physiology*condition*cladeThe random effect of species will be taken into account later.
# --------------------
# design cateogories
# --------------------
species<-as.character(unlist(design[1,]))
physiology<-as.character(unlist(design[2,]))
clade<-as.character(unlist(design[3,]))
condition<-as.character(unlist(design[5,]))
condition_physiology<-as.vector(paste(condition,physiology,sep="."))
condition_physiology_clade <- as.vector(paste(condition_physiology,clade,sep="."))
condition_physiology_clade <- as.vector(paste("group",condition_physiology_clade,sep=""))
cols<-colnames(counts)
ExpDesign <- data.frame(row.names=cols,
condition=condition,
physiology = physiology,
clade = clade,
species = species,
sample=cols)
ExpDesigndesign = model.matrix( ~0 + physiology*condition*clade, ExpDesign)
colnames(design)## [1] "physiologyFW"
## [2] "physiologyM"
## [3] "condition15_ppt"
## [4] "cladeClade2"
## [5] "cladeClade3"
## [6] "physiologyM:condition15_ppt"
## [7] "physiologyM:cladeClade2"
## [8] "physiologyM:cladeClade3"
## [9] "condition15_ppt:cladeClade2"
## [10] "condition15_ppt:cladeClade3"
## [11] "physiologyM:condition15_ppt:cladeClade2"
## [12] "physiologyM:condition15_ppt:cladeClade3"# check rank of matrix
#Matrix::rankMatrix( design )
#dim(design)3 Filtering and Normalization
Genes with low expression across samples were dropped from the analysis using a conservative approach. The function filterByExpr was used on the raw counts matrix. For each condition_physiology group (regardless of species), each sample must have a minium count of 10, and a group minium total count of 100. This reduced the counts table to the following dimensions (genes x samples):
gene.names<-rownames(counts)
counts<-as.matrix(as.data.frame(sapply(counts, as.numeric)))
rownames(counts)<-gene.names
#class(counts)
keep<-filterByExpr(counts,design = design,group=condition_physiology,min.count = 10, min.total.count = 100)
counts.filt <- counts[keep,]
dim(counts.filt)## [1] 21401 81#write.csv(counts.filt,"../../../21k_counts_filt_30April2019.csv")4 Genes of Interest
After filtration, we will check the counts matrix for the presence of several genes of interest. These genes have demonstrated responsiveness to salinity change from past studies.
| gene | Funhe | Ensembl | |
|---|---|---|---|
| zymogen granule membrane protein 16 | Funhe2EKm029929 | ENSFHEP00000007220.1 | |
| zymogen granule membrane protein 16 | Funhe2EKm029931 | ENSFHEP00000025841 | |
| solute carrier family 12 member 3-like (removed) | Funhe2EKm006896 | ENSFHEP00000009214 | |
| chloride channel, voltage-sensitive 2 (clcn2), transcript variant X2 (removed) | Funhe2EKm024148 | ENSFHEP00000019510 | |
| ATP-sensitive inward rectifier potassium channel 1 | Funhe2EKm001965 | ENSFHEP00000015383 | |
| inward rectifier potassium channel 2 | Funhe2EKm023780 | ENSFHEP00000009753 | |
| aquaporin-3 | ENSFHEP00000006725 | ||
| cftr | Funhe2EKm024501 | ENSFHEP00000008393 | |
| polyamine-modulated factor 1-like | Funhe2EKm031049 | ENSFHEP00000013324 | |
| sodium/potassium/calcium exchanger 2 | Funhe2EKm025497 | ENSFHEP00000034177 | |
| septin-2B isoform X2 | ENSFHEP00000015765 | ||
| CLOCK-interacting pacemaker-like | Funhe2EKm026846 | ENSFHEP00000017303 | |
| vasopressin V2 receptor-like | Funhe2EKm026721 | ENSFHEP00000000036 | |
| sodium/potassium-transporting ATPase subunit beta-1-interacting protein 1 | Funhe2EKm001174 | ENSFHEP00000031108 | |
| septin-2 | Funhe2EKm012182 | ENSFHEP00000016853 | |
| otopetrin-2 | Funhe2EKm035427 | ENSFHEP00000026411 | |
| claudin 8 | Funhe2EKm037718 | ENSFHEP00000006282 | |
| claudin 4 | Funhe2EKm007149 | ENSFHEP00000003908 |
If the Ensembl ID displays below, the gene is present in the whole data set and has not been filtered.
# ---------------------------
# Andrew Whitehead's genes of interest:
# ---------------------------
# zymogen granule membrane protein 16
# Funhe2EKm029929
# ENSFHEP00000007220.1
countsfilt <- counts.filt
countsfilt$row <- rownames(countsfilt)
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000007220.1"]
goi## [1] "ENSFHEP00000007220.1"# zymogen granule membrane protein 16
# Funhe2EKm029931
# ENSFHEP00000025841
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000025841"]
goi## [1] "ENSFHEP00000025841"# solute carrier family 12 member 3-like (removed)
# Funhe2EKm006896
# ENSFHEP00000009214
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000009214"]
goi## [1] "ENSFHEP00000009214"# chloride channel, voltage-sensitive 2 (clcn2), transcript variant X2 (removed)
# Funhe2EKm024148
# ENSFHEP00000019510
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000019510"]
goi## [1] "ENSFHEP00000019510"# ATP-sensitive inward rectifier potassium channel 1
# Funhe2EKm001965
# ENSFHEP00000015383
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000015383"]
goi## [1] "ENSFHEP00000015383"# inward rectifier potassium channel 2
# Funhe2EKm023780
# ENSFHEP00000009753
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000009753"]
# --------------------------------
# other salinity genes of interest
# --------------------------------
# ============================================
# aquaporin-3
# ENSFHEP00000006725
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000006725"]
goi## [1] "ENSFHEP00000006725"# cftr
# Funhe2EKm024501
# ENSFHEP00000008393
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000008393"]
goi## [1] "ENSFHEP00000008393"# polyamine-modulated factor 1-like
# Funhe2EKm031049
# ENSFHEP00000013324
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000013324"]
goi## [1] "ENSFHEP00000013324"# sodium/potassium/calcium exchanger 2
# ENSFHEP00000034177
# Funhe2EKm025497
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000034177"]
goi## character(0)# septin-2B isoform X2
# ENSFHEP00000015765
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000015765"]
goi## [1] "ENSFHEP00000015765"# CLOCK-interacting pacemaker-like
# ENSFHEP00000017303
# Funhe2EKm026846
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000017303"]
goi## [1] "ENSFHEP00000017303"# vasopressin V2 receptor-like
# Funhe2EKm026721
# ENSFHEP00000000036
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000000036"]
goi## [1] "ENSFHEP00000000036"# sodium/potassium-transporting ATPase subunit beta-1-interacting protein 1
# ENSFHEP00000031108
# Funhe2EKm001174
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000031108"]
goi## [1] "ENSFHEP00000031108"# septin-2
# Funhe2EKm012182
# ENSFHEP00000016853
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000016853"]
goi## [1] "ENSFHEP00000016853"# otopetrin-2
# Funhe2EKm035427
# ENSFHEP00000026411
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000026411"]
goi## [1] "ENSFHEP00000026411"# claudin 8
# Funhe2EKm037718
# ENSFHEP00000006282
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000006282"]
goi## [1] "ENSFHEP00000006282"# claudin 4
# ENSFHEP00000003908
# Funhe2EKm007149
goi <- countsfilt$row[countsfilt$row == "ENSFHEP00000003908"]
goi## [1] "ENSFHEP00000003908"all_goi<-c("ENSFHEP00000007220.1","ENSFHEP00000025841","ENSFHEP00000019510",
"ENSFHEP00000015383","ENSFHEP00000009753","ENSFHEP00000006725","ENSFHEP00000008393",
"ENSFHEP00000013324","ENSFHEP00000001609","ENSFHEP00000013324","ENSFHEP00000034177",
"ENSFHEP00000015765","ENSFHEP00000017303","ENSFHEP00000000036","ENSFHEP00000031108",
"ENSFHEP00000016853","ENSFHEP00000003908")5 Exploratory Plots
Log counts before normalization:
# log counts before DE
boxplot(log(counts.filt+1), las = 2, main = "")
Log counts after cpm normalization
# ---------------
# DE analysis
# ---------------
genes = DGEList(count = counts.filt, group = condition_physiology_clade)
genes = calcNormFactors( genes )
# write normalized counts
dir <- "~/Documents/UCDavis/Whitehead/"
tmp <- as.data.frame(cpm(genes))
tmp$Ensembl <- rownames(tmp)
tmp <- dplyr::select(tmp, Ensembl, everything())
write.csv(tmp, file = file.path(dir, "normalized_counts.csv"), quote = F, row.names = F)
vobj = voom( genes, design, plot=TRUE)
lcpm <- cpm(genes$counts, log = TRUE)
# log counts after DE
boxplot(lcpm, las = 2, main = "")
plot(colSums(t(lcpm)))
Voom weights:
vwts <- voomWithQualityWeights(genes, design=design, normalization="quantile", plot=TRUE)
The random effects of species are taken into account with the duplicateCorrelation function, which estimates the correlation between species. This reflects the between-species variability. The higher the correlation (0-1.0), the higher the variability between species.
corfit <- duplicateCorrelation(vobj,design,block=ExpDesign$species)
corfit$consensus## [1] 0.7596421#[1] 0.7513815.0.1 Individuals clustered by overall expression
counts_round<-round(counts.filt,digits=0)
dds <- DESeqDataSetFromMatrix(countData = counts_round,colData = ExpDesign,design = design)## converting counts to integer moderld <- vst(dds, blind = FALSE,fitType='local')
sampleDists <- dist(t(assay(rld)))
df <- as.data.frame(colData(dds)[,c("physiology","condition","clade")])
sampleDistMatrix <- as.matrix( sampleDists )
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors, annotation = df, show_rownames=F)
5.0.2 Individuals clustered by Top 100 genes
5.0.3 Individuals clustered by top 50 gene expression
5.0.4 PCA for overall expression
cowplot::plot_grid( plotPCA(rld, intgroup="condition"),
plotPCA(rld, intgroup="physiology"),
plotPCA(rld, intgroup="clade"),
plotPCA(rld, intgroup=c("clade","physiology","condition")),
align="c", ncol=2)
6 Contrasts
After running the lmFit function, which fits the linear model for each gene in the matrix and takes the random effects correlation into account, the resulting linear model fit is then used to compute the estimated coefficients and standard errors for a given set of contrasts:
fitRan <- lmFit(vobj,design,block=ExpDesign$species,correlation=corfit$consensus)
#colnames(coef(fitRan))
y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade,
data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition|clade*physiology ), reverse = T)
contrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
contrasts <- paste0(contrasts, "_", tmp$clade, "_", tmp$physiology)
rownames(contrast.matrix) <- contrasts
contrasts## [1] "15_ppt_v_0.2_ppt_Clade1_FW" "15_ppt_v_0.2_ppt_Clade2_FW"
## [3] "15_ppt_v_0.2_ppt_Clade3_FW" "15_ppt_v_0.2_ppt_Clade1_M"
## [5] "15_ppt_v_0.2_ppt_Clade2_M" "15_ppt_v_0.2_ppt_Clade3_M"tables <- lapply(contrasts, function(contr){
#print(contr)
cm <- contrast.matrix[contr,]
ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
tmp <- contrasts.fit(fitRan, contrasts = cm)
tmp <- eBayes(tmp)
tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
#tmp3 <- tmp2
#tmp3$row <- rownames(tmp3)
#tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
#tmp3 <- tmp3[order(tmp3$adj.P.Val),]
filename <- paste0(contr, ".csv")
#write.csv(tmp2, file = file.path(dir, filename),quote = F)
tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
header1 <- 6
names(header1) <- paste0("Top 20 genes for ", contr)
header2 <- 6
names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
header3 <- 6
names(header3) <- paste0("Output file is ", filename)
tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
tab <- tab %>% kable_styling()
return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}
)7 Three-way contrasts
The number of genes significant for the three-way interaction of condition:physiology:clade:
sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")| Comparison | Number of genes with adjusted P < 0.05 |
|---|---|
| 15_ppt_v_0.2_ppt_Clade1_FW | 13 |
| 15_ppt_v_0.2_ppt_Clade2_FW | 105 |
| 15_ppt_v_0.2_ppt_Clade3_FW | 46 |
| 15_ppt_v_0.2_ppt_Clade1_M | 533 |
| 15_ppt_v_0.2_ppt_Clade2_M | 272 |
| 15_ppt_v_0.2_ppt_Clade3_M | 105 |
8 Two-way contrasts
These genes demonstrate a conserved response to experimental condition across M or FW physiologies, regardless of clade.
y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade,
data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition|physiology ), reverse = T)## NOTE: Results may be misleading due to involvement in interactionscontrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
contrasts <- paste0(contrasts, "_", tmp$physiology)
rownames(contrast.matrix) <- contrasts
contrasts## [1] "15_ppt_v_0.2_ppt_FW" "15_ppt_v_0.2_ppt_M"tables <- lapply(contrasts, function(contr){
#print(contr)
cm <- contrast.matrix[contr,]
ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
tmp <- contrasts.fit(fitRan, contrasts = cm)
tmp <- eBayes(tmp)
tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
#tmp3 <- tmp2
#tmp3$row <- rownames(tmp3)
#tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
#tmp3 <- tmp3[order(tmp3$adj.P.Val),]
filename <- paste0(contr, ".csv")
#write.csv(tmp2, file = file.path(dir, filename),quote = F)
tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
header1 <- 6
names(header1) <- paste0("Top 20 genes for ", contr)
header2 <- 6
names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
header3 <- 6
names(header3) <- paste0("Output file is ", filename)
tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
tab <- tab %>% kable_styling()
return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}
)The number of genes significant for the two-way interaction of condition:physiology, independent of clade:
sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")| Comparison | Number of genes with adjusted P < 0.05 |
|---|---|
| 15_ppt_v_0.2_ppt_FW | 71 |
| 15_ppt_v_0.2_ppt_M | 616 |
9 Condition - main effect
These genes demonstrate a main effect of condition (15 ppt vs. 0.2 ppt), regardless of clade or physiology.
y <- rnorm(n = nrow(design))
dummy.mod <- lm(y ~ physiology*condition*clade,
data = ExpDesign)
pairs <- pairs(emmeans(dummy.mod, ~condition), reverse = T)## NOTE: Results may be misleading due to involvement in interactionscontrast.matrix <- pairs@linfct
tmp <- pairs@grid
contrasts <- gsub(" ", "", tmp$contrast)
contrasts <- gsub("-", "_v_", contrasts)
rownames(contrast.matrix) <- contrasts
contrasts## [1] "15_ppt_v_0.2_ppt"tables <- lapply(contrasts, function(contr){
#print(contr)
cm <- contrast.matrix[contr,]
ph <- sapply(strsplit(as.character(contr), "_"), tail, 1)
cl <- sapply(strsplit(as.character(contr), "_"), tail, 2)
tmp <- contrasts.fit(fitRan, contrasts = cm)
tmp <- eBayes(tmp)
tmp2 <- topTable(tmp, n = Inf, sort.by = "P")
#tmp3 <- tmp2
#tmp3$row <- rownames(tmp3)
#tmp3 <- merge(ann,tmp3,by.x = "ensembl_peptide_id", by.y = "row", all = TRUE)
#tmp3 <- tmp3[order(tmp3$adj.P.Val),]
filename <- paste0(contr, ".csv")
#write.csv(tmp2, file = file.path(dir, filename),quote = F)
tab <- kable(head(tmp2, 20), digits = 5, row.names = F)
header1 <- 6
names(header1) <- paste0("Top 20 genes for ", contr)
header2 <- 6
names(header2) <- paste0("Number of genes with adjusted P < 0.05 = ", length(which(tmp2$adj.P.Val < 0.05)))
header3 <- 6
names(header3) <- paste0("Output file is ", filename)
tab <- tab %>% add_header_above(header3, align = 'l') %>% add_header_above(header2, align = 'l') %>% add_header_above(header1, align = 'l', font_size = "larger", bold = T)
tab <- tab %>% kable_styling()
return(list(tab, nump = length(which(tmp2$adj.P.Val < 0.05))))
}
)The number of genes significant for the main effect condition:
sigps <- unlist(lapply(tables, function(x)x[[2]]))
sumtab <- data.frame(Comparison = contrasts, `Number of genes with adjusted P < 0.05` = sigps,
check.names = F)
kable(sumtab, align = 'c') %>% kable_styling() %>%
add_header_above(c("Overview of results" = 2), font_size = "larger", bold = T, align = "l")| Comparison | Number of genes with adjusted P < 0.05 |
|---|---|
| 15_ppt_v_0.2_ppt | 207 |
10 Heatmaps
10.1 Marine-Clade 1 (three-way) response
rld <- log2(mean_norm_counts_ordered_M_Clade1_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade1_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.2 Freshwater-Clade 1 (three-way) response
rld <- log2(mean_norm_counts_ordered_FW_Clade1_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade1_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.3 Marine-Clade 2 (three-way) response
rld <- log2(mean_norm_counts_ordered_M_Clade2_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade2_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.4 Freshwater-Clade 2 (three-way) response
rld <- log2(mean_norm_counts_ordered_FW_Clade2_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade2_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.5 Marine-Clade 3 (three-way) response
rld <- log2(mean_norm_counts_ordered_M_Clade3_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_Clade3_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.6 Freshwater-Clade 3 (three-way) response
rld <- log2(mean_norm_counts_ordered_FW_Clade3_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_Clade3_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.6.1 Marine physiology (two-way) conserved response
rld <- log2(mean_norm_counts_ordered_M_sig+1)
geneDists <- dist(mean_norm_counts_ordered_M_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.6.2 Freshwater physiology (two-way) conserved response
rld <- log2(mean_norm_counts_ordered_FW_sig+1)
geneDists <- dist(mean_norm_counts_ordered_FW_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
10.6.3 Condition response
rld <- log2(mean_norm_counts_ordered_condition_sig+1)
geneDists <- dist(mean_norm_counts_ordered_condition_sig)
df <- data.frame(ph,cl, condition,stringsAsFactors=FALSE)
rownames(df) <- colnames(rld)
pheatmap(rld, show_rownames=FALSE,
clustering_distance_rows = geneDists, cluster_cols= FALSE,annotation_col=df)
sessionInfo()## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] data.table_1.12.2 emmeans_1.3.4
## [3] kableExtra_1.1.0 biomaRt_2.38.0
## [5] vsn_3.50.0 lattice_0.20-38
## [7] gplots_3.0.1.1 edgeR_3.24.3
## [9] limma_3.38.3 DESeq2_1.22.2
## [11] SummarizedExperiment_1.12.0 DelayedArray_0.8.0
## [13] BiocParallel_1.16.6 matrixStats_0.54.0
## [15] Biobase_2.42.0 GenomicRanges_1.34.0
## [17] GenomeInfoDb_1.18.2 IRanges_2.16.0
## [19] S4Vectors_0.20.1 BiocGenerics_0.28.0
## [21] ggrepel_0.8.0 knitr_1.22
## [23] tidyr_0.8.3 dplyr_0.8.0.1
## [25] RColorBrewer_1.1-2 cowplot_0.9.4
## [27] pheatmap_1.0.12 gtools_3.8.1
## [29] ggplot2_3.1.1 MASS_7.3-51.4
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.4-1 estimability_1.3 htmlTable_1.13.1
## [4] XVector_0.22.0 base64enc_0.1-3 rstudioapi_0.10
## [7] affyio_1.52.0 bit64_0.9-7 AnnotationDbi_1.44.0
## [10] mvtnorm_1.0-10 xml2_1.2.0 splines_3.5.2
## [13] geneplotter_1.60.0 Formula_1.2-3 jsonlite_1.6
## [16] annotate_1.60.1 cluster_2.0.8 BiocManager_1.30.4
## [19] readr_1.3.1 compiler_3.5.2 httr_1.4.0
## [22] backports_1.1.4 assertthat_0.2.1 Matrix_1.2-17
## [25] lazyeval_0.2.2 acepack_1.4.1 htmltools_0.3.6
## [28] prettyunits_1.0.2 tools_3.5.2 coda_0.19-2
## [31] gtable_0.3.0 glue_1.3.1 GenomeInfoDbData_1.2.0
## [34] affy_1.60.0 Rcpp_1.0.1 gdata_2.18.0
## [37] preprocessCore_1.44.0 xfun_0.6 stringr_1.4.0
## [40] rvest_0.3.2 statmod_1.4.30 XML_3.98-1.19
## [43] zlibbioc_1.28.0 scales_1.0.0 hms_0.4.2
## [46] yaml_2.2.0 memoise_1.1.0 gridExtra_2.3
## [49] rpart_4.1-13 latticeExtra_0.6-28 stringi_1.4.3
## [52] RSQLite_2.1.1 highr_0.8 genefilter_1.64.0
## [55] checkmate_1.9.1 caTools_1.17.1.2 rlang_0.3.4
## [58] pkgconfig_2.0.2 bitops_1.0-6 evaluate_0.13
## [61] purrr_0.3.2 labeling_0.3 htmlwidgets_1.3
## [64] bit_1.1-14 tidyselect_0.2.5 plyr_1.8.4
## [67] magrittr_1.5 R6_2.4.0 Hmisc_4.2-0
## [70] DBI_1.0.0 pillar_1.3.1 foreign_0.8-71
## [73] withr_2.1.2 survival_2.44-1.1 RCurl_1.95-4.12
## [76] nnet_7.3-12 tibble_2.1.1 crayon_1.3.4
## [79] KernSmooth_2.23-15 rmarkdown_1.12 progress_1.2.0
## [82] locfit_1.5-9.1 grid_3.5.2 blob_1.1.1
## [85] digest_0.6.18 webshot_0.5.1 xtable_1.8-3
## [88] munsell_0.5.0 viridisLite_0.3.0